%4-sep-2007. This code reads Julia and ACe data sets
% pick days with cont. common measurements
% 

load c:\manoj\projects\ace\Julia_W_new w_climate Julia_W fday julia_euvac julia_solar_flux;
load c:\manoj\projects\ace\OMNI_ELEC_new;
load c:\manoj\geomag\indices\aplist.mat;

%%

w=w_climate; %This copies the Julia drift (measured)-Julia drift (modeled) to w
ace_fday = floor(ace_all(:,1));
ap_lower_limit = 0; %lower limt if Ap
solar_flux_limit = 0;
N_seg = 1;%the increasing counter for array JULI_SEG & ACE_SEG
N_data = 0;
phase_delay = 17;%minutes
mjd_date = datenum(2000,1,1);
plc = 'r';%plot color

phase_data = [0,10,17,25];

for i = 1:length(phase_data),
    legend_string(i,:) = sprintf('%4.1f',phase_data(i));
end;


for neeta = 1:length(phase_data),

clear ACE_SEG JULI_SEG TIME_SEG;
N_seg = 1;%the increasing counter for array JULI_SEG & ACE_SEG
N_data = 0;


ace_fday = floor(ace_all(:,1));

for i = 1: size(w,1),
   if sum(diff(Julia_W(i).k))+1 == length(Julia_W(i).k) &...
      sum(isnan(w(i,:))) <= 60 & julia_euvac(i) > solar_flux_limit, 

    L = ace_fday == Julia_W(i).fday;
    if sum(L) > 0,
        
    ace_time = ace_all(L,1)+(phase_data(neeta))/(60*24);%Try to advance 60 VS model/ ace time minutes
    ace_data = ace_all(L,2);

   
%
    L = isnan(ace_data);
    if sum(L) < 2 & sum(L) ~= 0, %NOTE more number of missing points, more phase shift at lower frequencies
        ace_data = interp1(ace_time(~L),ace_data(~L),ace_time);% in ace1, sum(L) <, 10 used, so corresponding phase shift to be expected
    end;%to make this phase difference similar to transfer function's phase difference, < < 2 should be used
    %bottom line use sum(L) < 10 (265 paisr)  for magnitudes , use sum(L) <
    %2 for phases (164 pairs)
%
   if abs(mean(diff(ace_time))-0.0035) <= 1e-004, %Use this with ACE min averages
        JULI = w(i,Julia_W(i).k);
        ACED = interp1(ace_time, ace_data, fday(i,Julia_W(i).k),'linear');
      

        L = fday_ap >= fday(i,Julia_W(i).k(1)) & fday_ap <= fday(i,Julia_W(i).k(end));
        mean_ap = mean(ap(L));

%         L = fday_gsm >= fday(i,Julia_W(i).k(1)) & fday_gsm <= fday(i,Julia_W(i).k(end));
%         mean_imf_bz = nanmean(ace_gsm_bz(L));
        

%       if sum(isnan(ACED)) < 1 & mean_imf_bz >=0,
         if sum(isnan(ACED)) < 1 & mean_ap >= ap_lower_limit,
%         if sum(isnan(ACED)) < 1 & mean_imf_bz <= 1 & mean_ap >=20,

            JULI_SEG(N_seg:N_seg+71) = JULI(1:72);
            ACE_SEG(N_seg:N_seg+71) = ACED(1:72);
            TIME_SEG(N_seg:N_seg+71) = fday(i,Julia_W(i).k(1:72));
            N_seg = N_seg+72;
            N_data = N_data+1;
            data_length(N_data) = length(JULI);
               
         end;
    end;
    end;
end;
end;

[Cxx,F] = mscohere(JULI_SEG,ACE_SEG,hanning(72),0,72,1/(5*60)); %1/(5*60) = sampling frequency 
[Pxy,F] = cpsd(ACE_SEG,JULI_SEG,hanning(72),0,72,1/(5*60));
phase = angle(Pxy);

phase_all_meenu(neeta,:) = phase;

end;
 phase_all_meenu1=unwrap(phase_all_meenu');
%correct=[-0.05;-0.15;-0.12;];
% correct=-0.15;
%phase_all_meenu1(2:4,:) = phase_all_meenu1(2:4,:) - repmat(correct,[1,neeta]);
phase_all_meenu1(2,2:3) = phase_all_meenu1(2,2:3) +0.15;


%  phase_all_meenu1(:,2) = phase_all_meenu1(:,2)+0.15;
% phase_all_meenu1(2,:) = phase_all_meenu1(2,:)+0.4784;

alpha = 0.95;
ci = 1- (1-alpha)^(1/(N_data-1));
L = Cxx(2:end) > ci;
%L(1) = 0; % make sure zero frequency is out
figure1=figure;
axes('Parent',figure1,'XTick',[0.1 1 10],...
    'XScale','log',...
    'XMinorTick','on');
set(gca,'FontSize',16);
box('on');
hold('all');
xlabel('period in hours');
ylabel('phase difference (degrees)');
hold on;
%semilogx((1./(3600*F(2:end))),phase_all_meenu1(2:end,:)*180/pi,'-','LineWidth',1);
semilogx((1./(3600*F(L))),phase_all_meenu1(L,:)*180/pi,'*-','LineWidth',2);

axis([0.1,10,-100,250]);
grid;

semilogx((1./(3600*F(L))),F(L)*(17*60)*2*180,'k','LineWidth',2);;
legend('No time delay','IEF delayed by 17 minutes','Linear phase shift for 17 minutes delay');